Vortex Lattice Melting in 2D Superconducting Networks and 

Films 



M. Franz and S. Teitel 

Department of Physics and Astronomy, University of Rochester, Rochester, NY 14627 

(February 1, 2008) 



Abstract 

i We carry out an extensive Monte Carlo study of phase transitions in 2D 

^ I superconducting networks, in an applied magnetic field, for square and hon- 

l> , 

! eycomb geometries. We consider both systems with a dilute vortex density 

O 

a^ 

O 

^\ , The dilute case gives the continuum limit as g — > oo, and serves as a model 



1/q, and dense systems near "full frustration" with vortex density 1/2 — 1/q. 



^ I for a uniform superconducting film. For this dilute case, we find a transition 

a, 

, temperature Tc ~ 1/q, at which the vortex lattice unpins from the network 

C , 

Q , and forms a "floating solid" phase. At a higher temperature T^, this floating 

O 

solid melts into a vortex liquid. We analyze the transition at according to 

•i-H . 

^ I the Kosterlitz-Thouless theory of dislocation mediated melting in 2D. While 

we find a discontinuous jump in the vortex shear modulus at which is 

consistent with this theory, we find (in opposition to this theory) that the 

transition is weakly first order, and we find no evidence for a hexatic liquid 

phase. For the case near full frustration, we find that the system can be 

described in terms of the density of defects in an otherwise fully frustrated 

vortex pattern. These dilute defects result in similar behavior as that found 

in the dilute vortex system, with pinned, fioating, and liquid defect phases. 
PACS number: 64.60-i, 74.50+r, 74.60-w, 74.76-w 
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I. INTRODUCTION 



Two dimensional (2D) periodic superconducting networks, and in particular arrays of 
Josephson junctions, have served as a convenient theoretical and experimental model system 
in terms of which one can study, in a well controlled way, the effects of thermal fluctuations 
and pinning, on vortex structures and phase coherence in 2D superconductorsB. Such 2D 
superconductors have received renewed attention recently with the observation that many 
of the high temperature superconductors consist of weakly coupled layers, and so for some 
range of parameters may display effectively two dimensional behavior!'!. One focus of this 
renewed interest has been concerned with the melting of the 2D vortex lattice, induced by 
an applied magnetic field, in a uniform continuous superconducting film. Controversy has 
resulted as to whether such a vortex lattice even exists at any finite temperature, or whether 
a vortex liquid is the only thermodynamically stable statei. In this work we address the 
thermodynamic behavior of vortex structures in 2D superconducting systems. Our focus will 
be on behavior in discrete periodic networks, however our results will also yield conclusions 
concerning the behavior of uniform films. 

Despite a decade of theoretical work, many fundamental questions remain unresolved 
concerning the nature of the phase transitions in 2D superconducting networks. When a 
uniform transverse magnetic field is applied, it induces a fixed density of vortices into the 
network, as in the mixed state of a type II superconductor. However, unlike a uniform 
superconductor, for which the ground state is a periodic triangular lattice of equally spaced 
vortices, the discrete network structure serves as an effective periodic pinning potential, 
which at low temperatures confines the vortices to sit at the centers of the unit cells of 
the network. This can result in novel vortex structures at low temperature, determined by 
the competition between the repulsive vortex-vortex interaction, and the periodic pinning 
potential induced by the networkiH^. Finding the ground state vortex structure for an 
arbitrary value of vortex density, for a given periodic network, remains an unsolved problem. 
The phase transitions at finite temperature have remained largely unexplored except for a 



2 



few of the simplest cases 

An early conjecture by Teitel and Jayaprakashi (TJ) argued that the superconducting 
transition in such networks would be governed by commensurability effects. If one measures 
the dimensionless vortex density / as the number of magnetic field induced vortices per unit 
cell of the network, they predicted that for rational / = p/q, the transition temperature 
would vary discontinuously as Tc{p/q) ~ While experimental evidence for high order 
commensurability effects has been reported in Josephson junction arrays,lii simulations by 
HalseyEl have challenged this conjecture for large q. A similar conjecture by TjH concerning 
the behavior of the ground state critical current, ic{f), has since been disproven in simula- 
tions by Lobb and co-workers and by Straley,0 who argue that as / varies, idf) has 
a lower non-zero limit determined by the single body effects of a non-interacting vortex in 
a periodic pinning potential; this conclusion has also been arrived at analytically by Val- 
lat and Beck0. However the validity of the TJ conjecture with respect to Tc(/), which is 
intrinsically determined by many body effects, has remained unresolved. 

In this paper we attempt to study the TJ conjecture systematically, by carrying out 
Monte Carlo simulations of superconducting networks for two special classes of vortex den- 
sity. We first consider the dilute case of vortex densities f = ^/q, q integer, for both square 
and honeycomb networks. This dilute case, as g — > oo, can equivalently be viewed as the 
continuum limit, in which the lattice spacing of the periodic network decreases to zero for a 
fixed areal density of vortices. Our results for this case therefore also address the problem of 
vortex lattice melting in a uniform continuous superconducting film. Secondly, we consider 
vortex densities / = l/2 — 1/g, close to full frustration, on the square network. 

Our results may be summarized as follows. For the dilute case with large q, the low 
temperature state is a Bravais lattice of vortices, with long range translational order, pinned 
commensurably to the periodic network. At a critical temperature Tc(l/g) ~ there is a 
sharp first order phase transition to a floating triangular vortex lattice, which is depinned 
from the periodic potential of the network. This floating lattice displays the algebraic 
translational correlations characteristic of a 2D vortex lattice in a uniform continuum. The 



depinning transition Tc(/) satisfies the TJ conjecture, and marks the loss of true d.c. super- 
conductivity in the network, due to the flux flow resistance which will result from drift of 
the unpinned vortex lattice. At a higher T^, which becomes independent of g as 1/g ^ 0, 
this floating vortex lattice melts into an isotropic vortex liquid. We analyze this transi- 
tion according to the theory of dislocation mediated melting in 2D, due to Kosterlitz and 
Thouless, Nelson and Halperin,lli and YoungS (KTNHY). While we find good agreement 
with certain predictions of this KTNHY theory, we find evidence that the second order 
melting transition predicted by KTNHY is pre-empted by a weak first order transition. 

For the close to fully frustrated case, / = 1/2 — the ground state is everywhere 
like that of / = 1/2 (a checkerboard pattern of vortices on alternating sites), except for 
a superimposed commensurate Bravais lattice of missing vortices, or "defects," so as to 
give the desired density / < 1/2. The transitions in this system are then governed by the 
behavior of these defects. Upon heating, there is first a depinning transition Tc(/) of the 
defect Bravais lattice into a floating triangular defect lattice; this depinning follows the TJ 
conjecture, Tc(/) ~ 1/g, and marks the loss of true d.c. superconductivity. At a higher 
Tm, the floating defect lattice melts into an isotropic defect liquid. Finally, at a higher T^/, 
there is an additional sharp transition representing the disordering of the vortices forming 
the / = 1/2 like background. 

The remainder of our paper is organized as follows. In Section II we present the theo- 
retical model used to describe the superconducting network, and its relation to a uniform 
superconducting film. We review the KTNHY theory of 2D melting, and discuss the ob- 
servables we measure and the methods we use to analyze our data. Finally we describe our 
Monte Carlo procedure. In Section HI we present our results for the dilute case / = 1/g on a 
honeycomb network. This corresponds to vortices on the dual triangular lattice of sites. We 
use finite size scaling to test in detail the predictions of KTNHY. In Section IV we present 
our results for the dilute case / = 1/g on a square network. In Section V we present our 
results for the dense case of/ = 5/llona square lattice, and infer the behavior for more 
general densities / = 1/2 — 1/g. In Section VI we present our discussion and conclusions. 
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II. MODEL AND METHODS OF ANALYSIS 



A. Model for a superconducting network 

A two dimensional superconducting network in a magnetic field, is described by the 
Hamiltonian, 

nm=J2Ui9,-9,-A,,) (1) 

where 6*, is the fluctuating phase of the superconducting wavefunction on node i of a periodic 
network of sites. The sum is over pairs of nearest neighbor sites, representing the bonds of 
the network, and 

A^ = {2n/<^o)J' A-dl (2) 

are fixed constants, giving the integral of the magnetic vector potential across bond (ij) 
($0 = hc/2e is the magnetic flux quantum). U{6) is the interaction potential between 
neighboring nodes, and its argument is just the gauge invariant phase difference across the 
bond. U{6) is periodic in 6 with period 27r, and has its minimum at 6* = 0. We will be 
interested here in the case of a uniform applied magnetic field V x A = B, transverse to the 
plane of the network. In this case, the sum of the Aij going counter clockwise around any 
unit cell of the network is constant, and determined by the magnetic fiux through the cell, 

J2Aj = 2nAB/^o = 2rrf, (3) 

cell 

where A is the area of a unit cell of the network. / therefore is the number of fiux quanta 
of applied magnetic field, per unit cell. 

For an array of Josephson junctions, the interaction potential in Eq.(|l]) is taken as 
U{6) = —Jocos{6). For a superconducting wire network, in the London approximation, a 
more appropriate interactioniil is given by the Villain function,ii defined by, 



oo 

^2 



-U{e)/T ^ J2 e--'o(e-27rm)V2T 



where we take ks = 1- 

For the Villain interaction, one can show by duality transformationJH that the Hamilto- 
nian of Eq. (|1|) can be mapped onto the following 2D classical Coulomb gas, 

^=^E(^^-/)nr.-r,)K-/), (5) 

where the sum is over all sites i, j of the dual lattice of the periodic network (i.e. the 
sites i in Eq.(|]) lie at the centers of the unit cells of the network), rii = 0,±1,±2,... are 
integer "charges" representing vortices in the phases 6i, and the magnetic field flux density 
is represented by the uniform background charge — /. V{r) is the lattice Coulomb potential 
in 2D, which solves the equation, 

AV(r) = -27r5r,o, (6) 

where is the discrete Laplacian for the network. For large separations, ^(r) ~ ln|r|. 
In mapping from the network Hamiltonian given by Eqs.(|I]) and (^, to the Coulomb gas 
Hamiltonian of Eq.(^, we have followed convention^ by rescaling the temperatures so that 
TcG = 7xy/27rJo, where Tcg refers to the temperature in the Coulomb gas model, and Txy 
refers to the temperature in the network (also referred to as a "uniformly frustrated" XY 
modeli). Henceforth, we will denote Tcg as simply T. 

Our simulations will be carried out in terms of this Coulomb gas problem, rather than 
in terms of the phases 6i. Although the Villain interaction may give quantitative differences 
when compared to the cosine interaction of a Josephson array, since the two functions have 
the same symmetry, we expect that they will display the same qualitative critical behavior .il 

For our simulations, we work with a finite LxL grid of sites, and apply periodic boundary 
conditions to the Laplace Eq.(|^) defining the Coulomb potential V{r). In this case, V can 
be explicitly calculated in terms of its Fourier transform.0 For a square network of lattice 
constant qq, one finds 

^^""^ " ^ ? 2 - cos(k ■ ai) - cos(k ■ as) ' ^''^ 



where = L^, {ai, a2}= {aox, ooy} ^^re the basis vectors, and the summation is over all wave 
vectors consistent with the periodic boundary conditions, i.e. the set {k} = {(mi/L)bi + 
{m2/L)h2}, with mi, 1712 = 0, 1, 2 ... L — 1, and with {bi, = {(27r/ao)x, (27r/ao)y} the 
basis vectors of the reciprocal lattice. 

For a honeycomb network, the charges nj sit on the dual triangular grid of sites, and the 
Coulomb potential is given by, 

^^^"^ ~ 2iV ^ 3 - cos(k • ai) - cos(k ■ aa) - cos(k ■ ag) ' 
where {ai, a2}={aox, ao(x/2 + v^y/2)} are the basis vectors, = a2 — ai, and the wave 
vectors are determined by {bi,b2} = {(27r/ao)(x — (l/-\/3)y), (27r/ao)(2/-\/3)y}. 

The k = terms in the summations of Eqs.(|^) and (§) will cause a divergence in V^(r). 
In real space, this is a reflection of the infinite self energy of a point charge. Configurations 
with infinite total energy will carry zero weight in the partition function sum, and may 
therefore be excluded. To keep the energy of the Coulomb gas finite, we therefore impose 
the condition of overall charge neutrality 

T.i^^ - /) = 0. (9) 

i 

If we define Nc as the total number of charges in the system, then Eq. (|^) gives 

N, = Y,ni = fN. (10) 

i 

Thus the density of magnetic flux quanta /, is equal to the density of charges (vortices) 
Nc/N. In the neutral system, the infinite self energies will exactly cancel, and in place of 
V{r) we can use only the nonsingular part of the Coulomb potential (^ and (||) defined 
by V'{r) = V{r) — V{r = 0). For a given system size, we evaluate V^'(r) by numerically 
performing the summations indicated in Eqs.(|^ and (^. 

The ground state will therefore be a periodic vortex structure consisting of A^^^ sites with 
Hi = +1 (all other sites having = 0), spaced as equally apart as allowed by the network 
geometry. Understanding the behavior of this vortex structure at finite temperature will be 
one of the main goals of this work. 



B. Relation to a uniform superconducting film 



The Coulomb gas model of the preceding section can also be used to describe the melting 
of the vortex lattice in a uniform continuous superconducting film. For a superconducting 
film, the states of the system can be described by a complex wavefunction, '?/'(r) = \ip(r)\e^^'^^\ 
As shown by Pearl, 13 for a of film of thickness d, provided the sample size is smaller than 
the transverse magnetic penetration length = X'^/d, the magnetic field will be essentially 
uniform and constant throughout the film. In this case, the states ?/'(r) will be weighted in 
the partition function sum according to the Ginzburg-Landau free energy. 



2e , , , 
-V A h/> 

2 C 



(11) 



with V X A = B a fixed constant. The mean field solution that minimizes T-Cgl['4^], is similar 
to that foundS in three dimensions: (z) there is a triangular lattice of equally spaced vortices 
in the phase ^^(r); {ii) the areal density of vortices is B/^o, with an average separation of 



O'v ~ y^o/B; {in) the size of the normal core of a vortex is determined by ~ 
where a = determines the B = mean field transition temperature; [iv) the mean field 
phase transition at finite B occurs when C,o ^ a^. 

To include fluctuations, one should now sum the partition function over all fluctuations 
of ip{r) about the mean field solution. In doing so, one common approach has been to 
make the London approximation. Here one assumes that, outside of the normal vortex 
core, the amplitude of the superconducting wavefunction is kept constant, and only the 
phase fluctuates, i.e. ip{r) = \ip\e^^^^\ The London approximation is expected to be good 
whenever the bare vortex core radius is very much smaller than the average separation 
between vortices, ^ ciy] by (iv) above, this corresponds to temperatures well below the 
mean field phase transition. 

Substituting ^{r) = \tp\e''^^^^ into Eq. ([Tl|) results, within additive constants, in the sim- 
plified free energy 
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(12) 



where Jo = ^q/IGit^Xx, and the integral is imphcitly cut off at the vortex cores. Eq.(|T2D 
is just a continuum version of the network Hamihonian, Eq.(0). Following Halperin and 
Nelsoni^ who considered the B = case, and Huberman and Doniach0 and Fisheri^ who 
considered the finite B case, we note that Eq. ( [T^ ) can be mapped onto a continuum Coulomb 
gas of logarithmically interacting charges. For finite B, this can be written@ in the form of 
a one component plasma on a uniform background charge density B/^o, 

n=^J c/Vc/V[n(r) - B/<^o]V{r - r')[n(r ) - fi/<l>o]. (13) 

Here n{r) = (l/27r)z • V x V6' is the vorticity in the phase of the superconducting wavefunc- 
tion, determined by singular integer vortices rii at positions r^, n{r) = J^i^^i^i''^ ^ ^i)- ^(i") 
solves the 2D Laplace equation, V^V = — 27r(5(r). 

The Coulomb gas of Eq.(^, introduced in the preceding section as a description for 
a network, can now be viewed as a discrete approximation to the continuum problem of 
Eq. (p!3D . For a fixed areal density of vortices, B/^q, we recover the continuum Eq.(|l^) from 
the discrete Eq.(^ as we take the network lattice constant Oq —* 0. Since the number of 
vortices per unit cell in the network is / ~ a^B/^Q, we see that the continuum is equivalent 
to the / — > limit. Thus by studying the melting of dilute vortex lattices in a network, we 
can also learn about the melting of a vortex lattice in a uniform superconducting film. As 
in the previous section, the mapping between the Coulomb gas and the superconductor is 
obtained by measuring the Coulomb gas temperature Tea in units of 2njQ = ^q/Stt^Aj,, i.e. 

Finally, we note that the melting of the 2D vortex lattice, described by the continuum 
Coulomb gas Hamiltonian of Eq.(0), has been treated within the general 2D melting theory 
of KTNHY. Within this theory, Fisheiiil has estimated that the melting transition occurs 
more than an order of magnitude below the mean field transition. This observation completes 
the self consistency of the argument for using the London approximation. 
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C. Review of the theory of 2D melting 



The analysis of our results will be guided by the ideas of the theory of defect mediated 
melting in 2D, developed by KTNHY. Although our results are, in many aspects, in 
opposition to this KTNHY theory, it still represents a useful starting point in exploring the 
phenomenon of 2D melting. 

For the 2D harmonic crystal on a smooth substrate (i.e. in the absence of any one-body 
potential) it is well known that fluctuations in the long-wavelength phonon modes, lead to 
a logarithmic divergence in the displacements of the particles, destroying translational long 
range order at any finite temperature. This is a consequence of the rigorous Mermin- Wagner 
theorem^! concerning long range order in 2D. The standard theory of elasticity shows how- 
ever, that despite the absence of translational long range order, the low temperature phase 
of such a crystal is characterized by a slow power-law decay of translational correlationsJ^i@ 
very different from the fast exponential decay that one would expect in the liquid. This 
phenomenon has been termed "quasi-long range" order, and we shall refer to such a phase 
as a "2D solid" . Based on the ideas of Kosterlitz and Thouless,lii that the melting of such 
a 2D solid would be nucleated by the unbinding of topological lattice defects. Nelson and 
Halperinlii and Youngil formulated a theory (KTNHY) which predicted that 2D melting 
would occur via two separate second order KT-like transitions. In particular, they predicted 
that the 2D solid with algebraic translational correlations would become unstable to the un- 
binding of dislocation defect pairs at a temperature Tm, and melt into a new phase called the 
hexatic liquid. This hexatic phase would be characterized by short range translational order, 
but quasi-long range six-fold orientational order. As the temperature is increased, KTNHY 
predicted that this quasi-long range orientational order would eventually be destroyed by 
the unbinding of disclination defect pairs, and at Tj > Tm, the hexatic liquid would melt 
into a normal (isotropic) liquid with short range orientational order. We summarize this 
scenario by writing down the long-range limiting behavior predicted for the translational 
and orientational correlation functions. 



10 



For translational correlations, 



1 in perfect crystal (T = 0) 



r 



in 2D solid (0 < T < T^) (14) 



e in hexatic or normal liquid (T > T^) 

where r^j = |rj — rj| is the separation between particles i and j, and G is a reciprocal lattice 
vector of the perfect (triangular) crystal at T = 0. TjG,{T) is a temperature dependent 
exponent, which for the 2D harmonic crystal can be expressed in terms its shear modulus /x 
and bulk modulus A, as 

WG£(3^ 

4vrM2/i + A) ■ ^^^^ 

In the 2D vortex lattice, the bulk modulus A is infinite because of the long-ranged nature of 
the Coulomb interaction. The expression for tiq,{T) thus simplifies to, 

^c(T) = (16) 

A key prediction of the KTNHY theory is that if Gi is a shortest reciprocal lattice vector, 
then r^Qj (T) takes a discontinuous jump qX Tfji to zero from the universal value, 

VG\iT,-) = 3. (17) 

In what follows, we will directly test this prediction. We wish to stress however, that the 
behavior of translational correlations in the 2D solid phase, as given by Eq.(|l^, is a general 
result of continuum elastic theory, independent of all assumptions concerning the mechanism 
of the melting transition. It is only the universal jump in ?7Q^(Tm), and the existence of the 
hexatic phase, which are specific predictions of KTNHY. 

The six-fold orientational correlation function, according to KTNHY, behaves as 



ae-^»i/«6 + in 2D sohd (0 < T < T„ 



r, 



-m{T) 



in hexatic liquid (Tm < T < Ti) (18) 



in normal liquid (T > Ti 
11 



where '(9(rj) is the angle of the bond from particle i to its neatest neighbor, relative to some 
fixed reference direction, a is a proportionality constant of order one, and gives the 
value of the long range orientational order expected in the 2D solid phase. The exponent 
r]Q^{T), describing the quasi- long range order of the hexatic phase, is predicted to have a 
universal jump to zero at Tj from the value rj^^iTf) = 4. 

In the 2D solid, a relation between ip'^ and the vortex shear modulus /i can be derived 
from continuum elastic theory,il 



0^ ~ exp 











= exp 


[ 2|G|2j 



(19) 



where A ~ 27i/ay is an ultraviolet cutoff [a^ is the average separation between particles). 
Since we will independently measure y?^ and r^G in our simulation, we will use this relation 
as a check of the consistency of our results. 

For a periodic superconducting network, we have discussed how the discrete substrate 
of the network serves to induce a periodic pinning potential for the magnetic field induced 
vortices. To treat this case, we are therefore interested in how the above 2D melting scenario 
is altered by the presence of a periodic substrate. We shall be interested in the situation 
where the period of the substrate is sufficiently small compared to the spacing between 
particles, so that the essential features of the defect mediated melting theory remain intact. 
This problem has been treated by Nelson and HalperinEl. The main result of such a "fine- 
mesh" periodic perturbation is the appearance of a new phase at low temperatures, in 
which the 2D solid is commensurably pinned to the substrate. This phase has true long 
range translational order, and we shall refer to it as the "pinned solid". At a certain 
depinning temperature Tc < Tm, there is a transition to a 2D "fioating solid" phase, where 
the solid decouples from the substrate, and translational correlations behave identically to 
those of a 2D solid on a uniform substrate; this triangular fioating solid may in general 
be incommensurate with the periodic substrate. Increasing temperature, the fioating solid 
is expected to melt at T^, via the dislocation unbinding mechanism, into a liquid phase. 
On a triangular substrate, this liquid will have a small (but finite) long-ranged six-fold 
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orientational order induced by the substrate, at all temperatures. There should, however, 
be a temperature Tj where (pf{T) shows a significant drop, reminiscent of the disclination 
unbinding transition on the smooth substrate. This drop should become increasingly sharper 
as the ratio of substrate period to particle separation becomes smaller. On the square 
substrate. Nelson and Halperin predict that there will be a sharp Ising transition at a 
Ti > Tm, where quasi-long range six- fold orientational order in the liquid vanishes, and only 
the long range four-fold orientational order induced by the substrate remains. This Ising 
transition can be viewed as a "ghost" of the hexatic to normal liquid transition, which would 
occur in the absence of the periodic substrate. The four-fold orientational order, induced by 
the substrate, again persists at all higher temperatures. 

To conclude, we note again that the properties of the 2D floating solid, described above, 
follow solely from continuum elastic theory, independent any particular theory of melting. 
It is the existence of the hexatic liquid phase that is a specific prediction of the KTNHY 
melting theory. However, as pointed out by Nelson and Halperin,^! it is always possible 
that a "premature" unbinding of disclination pairs may lead to a direct melting of the 2D 
solid into the normal liquid. Such a transition is then expected to be first order. In this 
case, the KTNHY prediction, Eq. (P^D , for the universal value of ?7q^ at melting, becomes a 
lower bound, i](2\(T^) > 3. Results from various numerical simulations and experimentsS 
indicate that this first order behavior might indeed be prevailing in the various 2D systems 
studied so far. 

D. Observables and finite size scaling 

We now show how the predictions of the preceding section translate into the behavior 
of observables which can be directly measured in our MC simulation. There are two key 
issues that we wish to investigate in the superconducting networks: (i) the transition from 
the superconducting to the normal state, and (ii) the melting of the magnetic field induced 
vortex lattice. For {ii), our goal is to test the KTNHY theory of 2D melting, and so we will 
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be interested in studying both the translational and the orientational order of the vortex 
lattice. 

The superconducting to normal transition, marked by the loss of superconducting phase 
coherence, is measured by the vanishing of the helicity modulus, T(T), which measures the 
response of the system to applying a net twist, or phase gradient, to the phases 9i in the 
Hamiltonian, Eq.([l|). For the Villain interaction of Eq.(H), the helicity modulus can be 
shownS to be identical to the inverse dielectric function of the corresponding Coulomb gas 
of Eq.(ll), T/Jo = e~^, where is defined in the usual way. 



Here = J2i exp(— ik- rj) is the Fourier-transformed charge density. The vanishing of 
signals an insulator to metal transition in the Coulomb gas. The free charges characteristic of 
the conducting phase correspond to freely diffusing vortices in the superconducting network, 
which are responsible for the loss of phase coherence.0 In the simulation, the k — > limit 
is approximated by averaging over the smallest allowed nonzero wave vectors. 

Information on the translational order in the vortex lattice can be extracted from the 
structure function 



which we evaluate for all allowed wave vectors k = (mi/L)bi + (m2/-L)b2 in the first Brillouin 
zone (BZ) of the reciprocal lattice to the real space dual lattice of the superconducting 
network. A 2D intensity plot of 5'(k) serves as a simple tool for visualization of the different 
phases in the system. In analogy to the conventional X-ray scattering images, we expect 
S'(k) to display a periodic array of sharp delta-function Bragg peaks in a state with long 
range translational order, and a set of smooth concentric rings in a normal liquid phase. 
A phase with quasi-long range translational order, characterized by algebraic translational 
correlations, will be distinguished by a regular array of algebraically diverging peaks of 
finite width.B A hexatic liquid phase should appear as a set of concentric rings with six-fold 
angular modulation. 




(20) 




(21) 
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Apart from providing the simple visualization described above, the scaling of the heights 
of the peaks in S'(k), as a function of system size L, will serve as a good quantitative indicator 
of translational correlations in the system. Combining the definition of 5'(k) in Eg. (pTj) with 
Eq.(p!^ (note that in Eg. (|2TD , rij = 1 on a site containing a vortex and = on a site 
without a vortex) one easily obtains 



S{G) 



1 in pinned solid (T < T^) (22a) 

L-vg{t) floating solid {T, < T < T^) (22b) (22) 

{^+/Ly in hexatic or normal liquid (T > T„i) (22c) 

The finite size scaling analysis of the translational order, that we present in Section III, will 
be based on the above relations. In particular, a comparison of Eqs.(p2l) with our MC data 
will allow us to extract the temperature dependent exponent tiq,(T) and test the KTNHY 
prediction regarding the universal jump in Tj^^{T,^). We shall also determine the correlation 
length i+{T) in the liquid phase. 

There is an independent way to extract the exponents r^c (and thus the vortex shear 
modulus /i) without having to use finite size scaling. One can instead, for a given system 
size, fit to the heights of the peaks 5'(G), as a function of |G|. For the low order peaks, 
this dependence is roughly Gaussian, as can be seen by combining Eq.(p2}3) with the ex- 
pression for r^G in Eq.(|T6]). For the higher order peaks however, we need to rederive this 
dependence, since the prefactor (which is not shown in Eq.(|2^)) becomes important. Sub- 
stituting Eq.(p!4|) into the definition of the structure function, Eq.(pll), and approximating 
the summations by integrations, we get 

5(G) = c / ciV(r/a„)-''« ~ 27rca;''° f^drr^-"^^, (23) 

J J ay 

where i? ~ L is a long distance cutoff, is the average separation between vortices, and c 
is a proportionality constant of the order one. The integral is easily evaluated 

27rc 



S(G)/L^: 



{R/a,)-''^ - {R/a, 



(24) 
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and c is determined by the requirement that S{0)/L'^ = 1. One can see that for r^G ^ 2 
formula (p2[b) remains vahd, but for r^c — 2 it breaks down. In practice, since r/Ci < 1/3 by 
the KTNHY bound, and any exponent r/G at a given fixed temperature can be written as 

VG = VG,-i\G\V\G,n (25) 

formula (P^b) will approximately hold for the three shortest reciprocal lattice vectors G 
only. To draw quantitative conclusions regarding the exponent rjc, one must use the more 
accurate relation of Eq . (^i]) . 

Information on the bond orientational order will be obtained by measuring the four-fold 
and six-fold orientational correlation 

y'p(T) = ^,T.(^''^''-''^), P = 4,6 (26) 

c ij ^ ' 

where the sum is over sites with non- vanishing charges Ui = +1, and -^i is the bond orien- 
tation angle defined in the previous section. For a finite system, one expects to see a sharp 
drop in fp{T) at the transition from an orientationaly ordered phase, to a disordered or 
possibly hexatic phase. One can deduce the scaling of the orientational correlation ipeiT) 
with system size L, by combining Eq.(PB[) with the KTNHY prediction of Eq.([T^) 

27ra(^^/L)2 + cp^ in pinned or fioating sohd (0 < T < T^) (27a) 
f6 ~ <, in hexatic liquid {Tm < T < Ti) (27b) (27) 

{S,(i/Lf in normal liquid (T > Ti) (27c) 

Relations (|27|a) and (p7|c) hold for L ^ (otherwise one must include corrections 
~ exp(— L/^e))- These scaling relations for ip^ will be used extensively in our analysis, 
to test for the existence of a hexatic phase in our model. 



E. Monte Carlo algorithm 



For the purpose of developing a fast MC algorithm, it is important to realize that the 
physical phenomena described in the previous section occur at temperatures which are about 
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one order of magnitude lowerEl than the ordinary Kosterhtz-Thouless (KT) transitioned in 
the zero magnetic field, / = 0, case. This imphes that the role of vortex-antivortex pair 
excitations is negligible in the temperature range that we study, and that in the simulation 
we can restrict ourselves to the excitations caused by movement of the vortices induced by 
the external magnetic field B. We have explicitly verified that the energy of an isolated 
vortex-antivortex pair (in the Coulomb gas language a pair of (+, — ) charges) is always 
Epair ^ kBT, and thus in practice such an excitation would never be accepted in the MC 
simulation. Consequently, our updating scheme is as follows. In each step, one charge is 
selected at random and moved to a different site within a radius ro, which is chosen so as 
to maximize the acceptance rate. We find that values tq ~ ay/2 are optimal. The energy 
change AE = Enew — Eoid is then computed, and the excitation is accepted or rejected 
according to the standard Metropolis algorithm: 

accept if e'^^^"^ > x, 

where x is a random number uniformly distributed on the interval [0, 1). Here and hence- 
forth, we work in units in which /c^ = 1. such attempts we will refer to as one MC 
sweep. At low temperature, we also made global moves, by attempting to shift entire rows 
of charges by one space. Such moves are meant to model long wavelength shear excitations, 
and help to accelerate equilibration near the vortex lattice melting transition. 

Due to the long-ranged nature of the Coulomb potential, the most time consuming 
operation is the evaluation of AE. From Eq.(|^) we find that the energy change for moving 
a charge from the site Ri to the site R2 is 

AE = -Y^ V'{Ri - r,)n, + ^ V'{K2 - r,)n, - V'{Ki - R2), (28) 

j j 

where we have used the fact that V'{—r) = V'{r) and V'{0) = 0. In this form, each 
evaluation of AE is a computation of the order Nc, as j sweeps through all the sites with 
nonvanishing charge rij = +1. To speed up this process, we use an algorithm developed by 
Crest .ii At each site of the lattice we define a potential due to all charges in the system 
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(29) 



j 



Now each evaluation of 



AE 



F(Ri) + F(R2)-\/'(Ri-R2) 



(30) 



requires computation of only 0(1). Naturally, each time the excitation is accepted, it is 
necessary to update F{ri) at all sites 



This is a computation of order N. However, since the acceptance rate in the interesting 
temperature range is very low (typically below 1%), this method is faster than the direct 
approach of Eq. (^) . 

Data is collected by heating the system up from the ground state. At each temperature we 
discard 30, 000 MC sweeps to equilibrate the system. Then, starting from this equilibrated 
configuration, we perform several (typically 4 — 6) independent runs of 100, 000 sweeps each 
to sample physical quantities. In some cases, when evaluating quantities at the temperatures 
close to the critical point, substantially longer runs are carried out. Errors are estimated from 
the standard deviation of these independent runs. To verify the consistency of our results, 
we also perform cooling from a random configuration at high temperature; no substantial 
hysteresis is found. 

All simulations were carried out on Sparc 10 workstations. The time needed to equilibrate 
the system and sample the physical quantities at a given temperature T was typically several 
hours (depending on size) using 100% of the single processor power. For example, it took 
approximately 3 hours to carry out 100,000 MC sweeps at a temperature close to melting, 
for a medium-sized system A^^ = 81 with density / = 1/49. Our longest run, to sample the 
energy distribution near the depinning transition, took 189 hours for 4 x 10^ MC sweeps on 
the largest system of Nc = 169 and / = 1/49. 



Fnewiri) = Foid{r,) - V'iri - Ri) + V (r^ - R2), ^ = 1, 2, . . . AT. 



(31) 
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III. SIMULATIONS ON THE TRIANGULAR GRID: DILUTE CASE 

A. The results 

In this section we report our results from simulations of the Coulomb gas Hamiltonian 
on a triangular grid of sites (corresponding to a honeycomb superconducting network), for 
the dilute hmit / ^ 1 (or equivalently Oq <^ a^). In this case, we expect that our discretized 
model will well approximate the continuum. Some of these results have been reported by 
us previously.^ The advantage of choosing a triangular grid is that, for a given system size 
L, one can always choose / in such a way so as to accommodate a perfect, commensurate, 
triangular vortex lattice in the ground state. By contrast, this is never possible on a square 
grid. It is convenient to choose / = with m integer, since then each system size of the 

form L = s -m, (s integer), will accommodate a triangular ground state with = fL"^ = 
vortices. We have studied systematically densities / = with m = 3 to 12, and fixed 

Nc ^ 100. The results of our investigation are summarized in Figure for sufficiently 
dilute systems (/ < 1/25) we find three distinct phases. At low temperatures the vortex 
lattice is in a "pinned solid" phase, locked to the underlying grid. Above a sharp depinning 
temperature Tc(/), the vortices are in a "fioating solid" phase, which then melts at Tm into 
a normal vortex liquid. The properties of these phases will be discussed below. For denser 
systems with / > 1/25, the two transitions at and merge, and there is only a single 
transition from a pinned solid into a liquid. 

For a simple visualization of the three phases, we show in Figure |^ intensity plots of 
S'(k) at various T, for the specific case of / = 1/49 and Nc = 63. We also display the 
amplitude of 5'(k) along the symmetry axis ky. For T = 0.003 (Figure ^), just below 
Tc{f), we see a regular array of 5-function Bragg peaks, indicating long ranged translational 
order induced by pinning to the triangular grid. The width of these peaks corresponds to 
the finite resolution of wave vectors allowed by our finite system. At T = 0.0065 (Figure 
Pd), just bellow Tm, we see a regular array of peaks, but they are now of finite width. We 
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will show that these peaks are consistent with the power law singularities characteristic of 
the algebraic translational correlations expected for a 2D floating solid phase. The heights 
of the peaks along the symmetry axis are well described by a Gaussian, as expected from 
Eqs. (p^ ) and (pS]). Thus, for Tc(/) < T < we do have a floating vortex lattice, as in 
the continuum limit. For T = 0.0075 (Figure ^), slightly above we see a rotationaly 
invariant structure, typical for a liquid with short range correlations. Thus for T > Tm, the 
floating vortex lattice has melted into a normal liquid. It is interesting to note that we see 
no sign of angular modulation in the rings above T^. One might expect such a modulation 
due to the long ranged six-fold orientational order induced in principle by the underlying 
triangular grid; if the grid was too fine for this effect to be significant, modulation might 
still be present if a hexatic phase existed just above T^. 

In Figure ^, we plot versus T the inverse dielectric function e~^(T), and the orientation 
order correlation (Pq{T), for / = 1/49 and Nc = 169 (one of the largest systems that we 
have studied). We see that e^^(T) vanishes at the depinning transition Tc(/), signaling the 
loss of superconducting phase coherence in the floating solid phase. This is just a reflection 
of the fact that an unpinned vortex lattice, in the presence of any applied d.c. current (no 
matter how small), will be free to drift transversely to the current, resulting in a finite linear 
"flux flow" resistance. Our results explicitly show that the absence of phase coherence in 
this k ^ sense, does not imply the absence of a well defined vortex lattice. Considering 
the orientation order, we see that ipe sharply drops at Tc{f), but remains finite up to the 
melting temperature T^, where it drops again sharply to nearly zero values. The smallness 
of ipe above indicates that the six-fold orientational long range order which is induced 
in principle by the triangular grid, is indeed a negligibly small effect at the densities we 
are concerned with. We shall discuss this point in more detail in the following section. 
In the Figure ^ we show the dependence of Tc(/) and on the vortex density /, as 
estimated from the behavior of e^^(T) and (Pq{T), and checked against the behavior of the 
structure function S'(k). We see that only for sufficiently dilute systems, / < 1/25, is 
there a floating solid phase; for / > 1/25 there is only a single transition from a pinned 
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solid to a liquid. As / decreases, Tc{f) vanishes linearly with /, consistent with the TJ 
conjectureil for the loss of superconducting coherence. however, quickly approaches a 
finite constant Tm = 0.0070 ± 0.0005. In terms of the superconductor temperature, this 
means a vortex lattice melting at = 0.0070 ^Q/8n'^X±. This is well within the bounds 
0.0046 <Tm< 0.0086 estimated by Fisheril from the KTNHY theory. 



B. The melting transition: finite size scaling analysis 

To investigate if the melting transition at is indeed consistent with the KTNHY 
theory, we have carried out a detailed finite size scaling analysis for the density / = 1/49. 
This density has been chosen for two reasons. First, the estimated Tm is well separated 
from the depinning temperature T^, and hence the floating solid phase exists in a relatively 
wide interval of temperatures. Second, the density is not too small, and thus we are able to 
study systems with as many as A^^^ = 169 vortices. More dilute systems, with comparable 
Nc would require sizes that are currently out of reach of the computer power available to 
us. We have carried out extensive simulations for the system sizes L = 28, 35, . . . 91, and we 
have analyzed the size- dependencies of various physical quantities at several temperatures 
below and above T^. Our results follows. 

In Figure || we plot S'(Gi)/L^ as a function of L, on a log-log scale, for several different 
temperatures. Data for each temperature fall on a straight line, confirming the expected 
power-law behavior of Eg . (^2]) . These straight lines fall into three distinct groups. For 
T <Tc — 0.0045, S{Gi)/L?' ~ 1, indicating the long range order of the pinned lattice. For 
T^<T <Tm^ 0.007, we find algebraic decay, S{Gi)/L'^ ~ For T > T^, we find 

S'(Gi)/L^ ~ L~^, with X — > 2 as T increases, consistent with the short range order of a 
liquid. 

Thus, our data for the floating solid phase are consistent with the predictions of the 
continuum elastic theory, given by Eq. (p^), and in particular we may fit our data to this 
expression to obtain the translational correlation exponent ?7gi • We show our results in Table 
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|. We can now make quantitative comparison with the KTNHY theory, by noting that ?7gi 
first exceeds the KTNHY universal value of 1/3 (see Eq.(|T^ at T = 0.0065, very close to 
the melting temperature Tm — 0.0070 as estimated from the behavior of the orientational 
correlation ipe{T) of Figure |^. The slopes of the lines in Figure ^ also show an apparent 
discontinuous jump at this same T^. 

As a consistency check, we have also computed rjc^, where G2 = 2Gi. Using similar fits 
to S'(G2) as in Figure |, we determine the exponent rjc^, and show the results in Table |. 
We see that rjG^ ~ 4?7gi as expected, since according to Eq. ([T9D rjG ~ |Gp. 

As an alternative way of calculating r/Ci, we fit to the heights of the peaks in 5'(k) at 
all the available G, for a fixed size system, as described in Eqs.(p3D — (pSj). We found that 
the results are only weakly dependent on the precise value of the cutoff R of Eq.(p^); we 
therefore take R = L, which results in an excellent fit. We show one example of such a fit 
in Figure The exponent riG,-^{T), obtained in this way, is shown in Table for the sizes 
L = 63, 77, 91, with / = 1/49. We note, that despite a certain tendency to overestimation, 
these exponents are in reasonable agreement with those obtained from the finite size scaling. 
This method of extracting the shear modulus should be useful in situations where a finite 
size scaling analysis is difficult to handle, such as in systems with a large unit cell in the 
ground state. 

Let us now consider the orientational order. In Figure |^ we plot the orientational cor- 
relation ife{T) as a function of L for several temperatures. In the pinned solid, T < T^, 
tfe{T) — s> 1 as L increases, confirming the expected long-ranged orientational order of the 
perfect pinned triangular lattice. More interestingly, ipe(T) also approaches a finite value 

in the floating lattice phase, Tc < T < Tm, in agreement with continuum elastic theory. 
The solid lines in Figure |^ are from least squares fits to Eq. (p7|a) . The resulting fitted values 
of ip'^ are shown in Table |. Above Tm we attempt to fit to the power law of Eq.(|27[b) for 
a hexatic liquid, but we always find that using Eq.(P7|c) for an isotropic liquid, results in 
a distinctly better fit. Since the underlying triangular grid will in principle result in long 
range six-fold order at all temperatures, we have also fit our data above Tm to the form of 
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Eg. (p7|a) , which differs from Eg. (p7|c) only in the constant <y9^. As shown in Table | however, 
we always find ~ 0. Thus the discrete grid is playing a negligible role in the orientation 
order. To compare our fits above T^, we note that in most cases the x^-parameter of the 
fit to Eg. (p7|a) is 5 to 10 times smaller than that of Eg. (p7|b). The former fit is also much 
more stable in the sense that fitted parameters do not change significantly when the data is 
restricted to different ranges of L. We may therefore conclude that, in agreement with our 
investigation of the structure function, the floating solid melts directly into a normal liguid. 
The hexatic phase is either absent in our system, or it occurs only in some extremely narrow 
interval of temperatures, which makes it difficult to detect by numerical simulation. 

In order to see this another way, in Figure ^ we plot versus temperature our values of 
(f^{T) and T/rjG-^iT), obtained from our finite size scaling analysis. From Eg.(|l^) we see 
that T/riG^{T) is just proportional to the vortex lattice shear modulus /i. We see, that 
V?^(r) starts to drop at the same temperature that T /riG^{T) first drops below the KTNHY 
universal value of 3T (see Eg. (p^T]) ) , i.e. the temperature at which the floating solid starts 
to melt. The temperature range over which decays to zero is identical to the range over 
which T /riG^{T) decays. This suggests that the small but finite values of ip'^ which we find 
above are just a finite size effect, rather than a signature of the hexatic phase. Let us 
also note that the exponent riG^{T) has a physical meaning only below T^. Above Tm it 
is strictly infinite and we use it here, with some abuse of notation, simply as the exponent 
resulting from the fit of our data to the power law form of Eg.(l22|b). 

Having obtained the values of rjG,^{T) and y^Jf (T), we can now test the relation. Eg . , 
between the orientational long range order and the vortex lattice shear modulus /i, that 
should hold in the floating solid phase. Expressing /x in terms of the exponent r^Ci, Eg.(p!9D 
gives 

^o{T) = -K\n[ipt{T)l (32) 

with K = 2|Gp/9A^, where A = X{27T/ay), and A is a dimensionless constant of order unity. 
In Figure I we plot VgHT) and -1 / K \n[<p^ (T)] versus temperature. For K = 0.33, which 
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corresponds to A = 0.94, the two data sets lie on top of each other for all T below Tm, 
providing yet another consistency check for our calculation. 

By fitting our data above to Eqs. (p2|c) and (p^c) we have also extracted the correlation 
lengths ^+(T), associated with translational order, and i^iT), associated with six- fold bond 
orientational order. We are able to determine these from finite size scaling only up to an 
overall multiplicative factor. We determine this factor by assuming that at high T, the 
correlation lengths are equal to the average spacing between vortices, i.e. ^{T oo) = a„. 
With this assumption, i+{T) and ^6(^) are displayed in Figure 0. We see that both 



correlation lengths rapidly increase around T ~ 0.007. It is also evident from Figure [T^, 
that orientational correlations persist out to longer distances, up to higher temperatures, 
than translational correlations. 

C. The order of melting transition 

The absence of the hexatic phase, as deduced from our analysis of the orientational 
correlations, suggests the possibility that the transition is not of the KTNHY type, but is 
due to some other mechanism, such as domain wall proliferation. It might also be, that 
the unbinding of disclinations occurs simultaneously with the unbinding of dislocations. 
Such a possibility has been suggested in Ref. In any case, it is useful to determine 
the order of this melting transition. To examine the possibility that the transition is first 
order, we have used the histogram method due to Lee and Kosterlitz0. For various system 
sizes at / = 1/49, we measure the energy distribution P{E) ~ Q-^i^)/'^ near the melting 



temperature T^. In Figure [TT] we plot the resulting free energy F{E) versus E. Although 



our data are somewhat noisy, we see a clear double well structure with an energy barrier 



AF between two coexisting phases. The inset to Figure |Tl| shows the dependence of AF on 
the system size L. The energy barrier AF grows with L, strongly suggesting a first order 
transition. Our system sizes remain too small to see clearly the predicted scaling AF ~ L. 
For all sizes, the data have been taken at T = 0.0065, and then the energy distribution 
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is extrapolated, using the method of Ferrenberg and SwendsenJ^ to that temperature which 
gives two minima of equal depth. This criterion gives an improved estimate of the melting 
temperature, = 0.0066. A total of 10^ MC sweeps were performed for each size to measure 
the energy distribution P{E). We have checked the consistency of these measurements by 
calculating the energy of the system at various temperatures (above and below T^) using 
the extrapolated distributions P{E,T). We then compared these with energies obtained by 
direct simulation at those temperatures, and found good agreement for all temperatures not 
too far from T^. 

To conclude, the histogram method provides strong evidence that the melting transition 
is first order. This is consistent with our observation that the 2D solid melts directly into 
an isotropic liquid. The transition is weakly first order however, as can be seen from our 
result that the jump in riG-^^(T~) (and hence the vortex lattice shear modulus /i) at melting 
remains very close to the KTNHY universal value. 



D. The depinning transition 

Finally, we consider the order of the depinning transition at Tc(/). In their work on 
2D melting on a periodic substrate. Nelson and Halperinlii studied this "commensurate 
to floating" transition using renormalization group techniques. They concluded that the 
transition is most likely second order, with properties very similar to the floating solid 
to liquid melting transition discussed in Section IIC. To test this prediction, we use the 
histogram method applied at the depinning transition Tc(/), just as we did in the preceding 
section for melting at T^. Measuring the energy distribution P{E) at Tc{f), for / = 1/49 
and various L, we show the free energy F{E) versus E in Figure [1^. As was seen at Tm, 
we now similarly see a pronounced double well structure with barrier AF growing with the 
size of the system (see the inset). Again, this is a clear indication that the transition is 
flrst order. Due to the low acceptance rates at these low temperatures, we had to perform 
as many as 4 x 10^ MC sweeps for each system size, in order to get reasonably accurate 



25 



energy distributions. By finding the temperature that produces minima of equal depth, we 
estimate that T^if = 1/49) = 0.0046. 

IV. SIMULATIONS ON THE SQUARE GRID: DILUTE CASE 

A. The ground state 

In the present section we shall investigate the dilute limit (/ <^ 1) of the Coulomb gas, on 
a square grid of sites. This corresponds to a square periodic superconducting network, which 
has been the predominant geometry in experimental and theoretical studies of networks.BiHi 
Qualitatively, we find similar behavior as found in Section III for the triangular grid: a 
depinning transition Tc(/), from a commensurate pinned solid to a floating solid, followed 
at higher temperature by a melting transition Tm to a liquid. 

While the case of a square grid is more relevant to the physics of superconducting arrays, 
it is somewhat more difficult to study theoretically than the triangular grid of Section III. The 
main reason for this is the rich variety of ground state configurations that one can encounter 
for various system sizes and vortex densities /. The most extensive enumeration of such 
ground states, for both dilute and dense /, has been carried out by Straley and Barnett.i This 
richness in ground state structure is due to the intrinsic competition between the repulsive 
vortex-vortex interaction, which prefers the formation of a perfect triangular lattice, and 
the geometrical constraints implied by the presence of the square grid. Since a triangular 
lattice is incommensurate with a square grid, for small / the resulting ground states form 
high order commensurate approximations to a triangular lattice, that vary substantially as / 
varies. Thus, while in the triangular grid a density / = 1/ can always fit commensurably 
in a lattice of size L = s ■ m {s integer), for the square grid, a density of / = 1/g (g integer) 
will require a lattice of at least L = s ■ q to contain the commensurate ground state. It thus 
becomes too difficult to carry out detailed finite size scaling calculations at small /, as the 
lattice sizes needed quickly become too large to simulate. We therefore must be content with 
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a more qualitative analysis based on simulations at a fixed size system. A second problem, 
related to the high order commensurability of the ground state, is the existence of excited 
states that are nearly degenerate in energy with the ground state. This can sometimes cause 
equilibration problems, or leave uncertainty as to the configuration of the true ground state. 
Fortunately, these difficulties occur only at low temperatures, below the depinning transition 
Tc(/), where commensurability effects are crucial. In sufficiently dilute systems, the melting 
of the floating solid phase at is largely unaffected by such difficulties, and we find results 
familiar to the preceding section. 

We have performed simulations for systems with a wide spectrum of densities f = l/q {q 
integer) with 10 < g < 90. From inspection of the inverse dielectric constant e~^(T) and the 
orientational correlation (Pq{T), for a fixed size L, we estimate the depinning and melting 
transition temperatures, Tc(/) and T^, and we plot these values versus / in Figure 0. We 
see that above / ~ 1/30 the depinning and melting transitions merge, and there is only a 
single transition from pinned solid to liquid. Due to the varying commensurability of the 
ground state as / varies, the values Tc{f) and no longer decrease monotonically with 
/, as was found for the triangular grid. Nevertheless, we see that Tc(/) still tends linearly 
to zero as / decreases (dashed line), in agreement with the TJ conjecture. appears to 
saturate around 0.007, in agreement with the melting temperature found for the triangular 
grid. 

In order to find the ground state of the system for a given density / and size L, we have 
devised a simple program that scans all possible periodic vortex configurations, consistent 
with periodic boundary conditions, and evaluates their energy. When translational and in- 
version symmetries are accounted for, the total number of distinct configurations is relatively 
small, and even for the largest of the systems that we considered it took only few minutes 
to execute the program. The lowest energy configuration obtained in this manner was then 
taken as a candidate for the ground state. In many cases we have verified that this indeed 
was a true ground state by performing a slow MC cooling from a random configuration at 
high temperature. In all cases we found that for / = 1/q, the ground state has a q x q peri- 
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odicity. In Figure we display two typical examples of these ground state configurations. 



The almost perfect triangular lattice (a/68 x v68 x v72) in Figure 0a is for / = 1/60. 
Figure 0b shows the example of a nearly square vortex lattice (v^50 x ^/53 x ^/89) with 
/ = 1/51. In what follows we shall concentrate on these two special cases as representatives 
of two classes of systems with slightly different physical properties. 

B. Systems with "nearly triangular" ground state: / = 1/60 

Not very surprisingly, systems with an almost triangular ground state, such as / = 1/60 
shown in Figure [l^, behave in a fashion similar to systems on the triangular grid studied 
in Section III. We display the behavior of the inverse dielectric function e^^(T) for / = 1/60 
and L = 60 in Figure |15|a. In Figure p!5| b we show the six-fold and four-fold orientational 
correlations, y:>6(T) and f^iT). A sharp drop in e^^(T) around Tc(/) ^ 0.0045 signals the loss 
of superconducting phase coherence. Above Tc(/), is zero, but v^6(^) stays finite. Based 
on our experience from the triangular grid, we take this as a signature of a floating triangular 
solid with long-range orientational order. Tc(/) is thus a transition from a commensurate 
pinned solid, to an incommensurate floating solid. Around ~ 0.0075 we see that feiT) 
drops again to very small values; we take this as a signal that the floating solid has melted 
into a vortex liquid. 

In order to confirm this scenario, we calculate the structure function 5'(k) at various 
temperatures, and display the resulting intensity plots in Figure |16|. We clearly see the 
pinned solid (Figure |l6|a), the floating solid (Figure 0b), and the liquid (Figure |l6|c) phases. 
It is interesting to note that the rotational symmetry of the pinned and floating solids break 
the four-fold rotational symmetry of the square grid, leading to two possible degenerate 
orientations. In the liquid however, we see that the four-fold symmetry of the square grid 
is restored, with a strong four-fold angular modulation of the circular intensity peaks. This 
observation is also confirmed by a direct measurement of v^4(T) (see Figure iMlb), which is 



close to zero in the floating solid phase, but then rises sharply at the melting transition and 
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only slowly vanishes with increasing temperature. The small values of iPa{T) for Tc{f) < 
T < Tm are an indication of the extent to which the commensurate, slightly distorted 
triangular lattice of the ground state, beomes an incommensurate perfect triangular lattice 
in the floating solid phase. 

Since we are unable to carry out finite size scaling, we are unable to search in detail 
for the hexatic phase, or for the predictedlii Ising transition from the hexatic to the normal 
liquid. However, as the four-fold symmetry appears to be restored at the same temperature 
as the melting transition, we suggest that, as was found for the triangular grid, the hexatic 
phase is absent and the melting transition is first order. 

Although finite size scaling is not possible, we can nevertheless still obtain the transla- 
tional correlation exponent ?7gi by analyzing the decay of the peaks in the structure function, 
using the method discussed in connection with Eqs.(p^)-(|25|). In the present case the imple- 
mentation of this method is somewhat trickier than it was for the triangular grid, since, due 
to the incommensurability of the floating solid, the peaks in ^(k) do not have well defined 
positions G on the square reciprocal lattice. We overcome this difficulty by numerically 
scanning S'(k) for local maxima at a given distance from the center of the reciprocal lattice, 
and averaging over the heights of peaks of the same order. The peak heights estimated in 
this way are shown in Figure |l^ for several temperatures T in the floating lattice phase. 



Dashed lines are least square flts to the formula (|24D , and the extracted exponents ric^iT) 



are summarized in Table |T|. The accuracy of the flt appears to be as good as in the case 
of the triangular grid, and we therefore have good reason to believe that our determination 
of rjG^{T) (and thus the vortex shear modulus /x) is reasonably accurate. Once again we see 
from Table |T| that rjc^iT) flrst exceeds the universal KT value of 1/3 at T ~ 0.007, very 
close to the melting temperature — 0.0075 estimated from the behavior of (^^{T). 
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C. Systems with "nearly square" ground state: / = 1/51 



We now briefly describe our results from simulations on a system with density / = 1/51, 
which possesses a nearly square ground state (Figure p!4|b). Our results are for a system of 
size L = 51. The temperature dependence of the inverse dielectric function e~^(T) is shown 
in Figure |T^a. From the data, we estimate the depinning temperature to be Tc(/) ^ 0.0035. 
We note that even though / = 1/51 here is larger than the / = 1/60 studied in the 
previous section, we find 0.0035 = Tc(l/51) < Tc(l/60) = 0.0045, thus illustrating the 
nonmonotonic behavior of Tc(/) for small /. The significantly lower depinning temperature 
in the present case may be qualitatively understood as a result of the larger distortion 
of the ground state away from the perfect triangular lattice favored by the vortex-vortex 
interaction. This large distortion, which is favored by the pinning energy, comes at a cost 
in vortex-vortex interaction energy. The result is a reduced free energy difference between 
the pinned "distorted triangular" solid and the perfect triangular floating solid, and hence 
a reduced depinning temperature. A similar observation holds for the other values of / we 
have studied: systems with relatively lower Tc(/) compared to other nearby values of /, are 
those with greater distortion of the ground state from a triangular lattice. 

In Figure |T8|b we plot the temperature dependencies of the four-fold and six-fold ori- 
entational correlations f^iT) and (peiT). In contrast to the previously considered cases, 
the melting transition is barely visible here: one sees only a small kink in <f4{T) and an 
inconspicuous dip in feiT) near T = 0.005. For a clearer picture of melting, we show 
the structure function S(k) in Figure We see again the pinned solid (Figure [l9|a), the 
floating solid (Figure |l9|b), and the liquid (Figure |l9|c). Note that the peaks in the floating 
solid occur at distinctly different wavevectors k than in the pinned solid; this emphasizes 
the fact that Tc{f) is truly a transition from a commensurate "nearly square" lattice, to 
an incommensurate floating triangular lattice. Inspection of these intensity plots gives a 



melting transition of Tm — 0.005, in agreement with the value hinted at in Figure |T8|b. 
We note that this value is significantly lower than the value of 0.007 found in other cases. 
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Thus commensurability effects can also significantly lower the melting temperature. Such 
commensurability effects presumably become less significant as / decreases and the ground 
state becomes increasingly closer to a triangular lattice. We thus expect that the melting 
temperature becomes — 0.007 in the asymptotic / limit, as is indeed suggested by 
Figure 113. 



V. SIMULATIONS ON THE SQUARE GRID: NEAR FULL FRUSTRATION 

The square superconducting network with / = 1/2 has been the focus of extensive 
theoretical study in recent yearsiSil. As the Hamiltonians, Eq.(|l]) and Eq.(^, are periodic 
in / with period 1, / = 1/2 represents the strongest magnetic field, and most dense vortex 
configuration, discernable by the network. Thus this case is usually refered to as "fully 
frustrated." The ground state of this configuration is in some sense the simplest of all 
/ > 0, consisting of a checkerboard pattern of vortices, with rii = 1 and rij = on the 
two alternating sublattices of the square grid. This dense vortex lattice melt^ directly 
into a vortex liquid at Tm(l/2) ~ 0.13. Superconducting coherence vanisheJ^l at Tc(l/2) ~ 
T^(l/2). 

We now wish to study how the system behaves as / is varied slightly away from 1/2, 
in order to test the discontinuous behavior predicted by the TJ conjecture. We study in 
particular systems with f = 1/2 — 1/q, with integer q large. While the ground states for 
densities of this form have been studied by Straley and co-workers,0@ finite temperature 
properties have remain unexplored. 

A. / = 5/ll 

We first consider the particular case of / = 5/11, which may be written as / = 1/2 — 1/22. 



The correct ground state for this case, which we show in Figure 20, was first found by 



Kolachi and Straley.0 It consists of a periodic superlattice of vortex vacancies superimposed 
on an otherwise uniform / = 1/2 like background, and is periodic with a 22 x 22 unit cell. 
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Our motivation is to see whether or not this superlattice of vacancies (or "defects") can 
melt independently of the / = 1/2 like background, and if so, whether the resulting liquid of 
vacancies destroys superconducting coherence. Our analysis is similar to that in the previous 
section. 

Heating from the ground state, we show sample intensity plots of the structure function 



S(k), at different temperatures, in Figure Figure |21| a shows the low temperature phase 
at T = 0.010. The bright Bragg peaks at k = (±7r/ao, ±7r/ao) originate from the vortex 
ordering in the / = 1/2 like background, while the periodic Bravais lattice of less intense 
Bragg peaks is due to the defect superlattice, which at this low temperature is pinned to the 
substrate. Thus we have a "pinned defect solid" phase. As the temperature is increased. 



we find that the defect superlattice melts at ~ 0.015. In Figure ^T]b we show the system 
at T = 0.018, just above this melting. The defects no longer give rise to Bragg peaks, but 
instead we see the circular rings (with strong four-fold angular modulation) characteristic of 
a defect liquid. However the bright Bragg peaks at k = (±7r/ao, ±7r/ao) remain, indicating 
that the / = 1/2 like vortex background remains ordered. Upon increasing the temperature 
further, the ordered / = 1/2 background is also eventually destroyed at T^/ ~ 0.040. In 
Figure E^c we show the system at T = 0.055, above T^/. The peaks at k = (±7r/ao, ±7r/ao) 



have broadened to finite width, indicating the disordering of the / = 1/2 like background. 



To see the melting transitions more clearly, in Figure ^ we plot versus temperature the 
peak heights 5'(q*) and S'(Gi), with q* = (vr/ao, vr/oo) giving the ordering of the / = 1/2 like 
background, and Gi the shortest reciprocal lattice vector of the defect superlattice. We see 
that ^(Gi) vanishes sharply at — 0.015, where the pinned defect superlattice melts into 
a defect liquid. 5'(q*), however, remains at its T = value of unity for all temperatures up to 
T ^ 0.020, clearly demonstrating that the / = 1/2 like vortex background remains ordered 
throughout the defect melting transition, ^(q*) starts to drop to zero around T^' — 0.04, 
where, based on the structure function intensity plots, we have estimated that the / = 1/2 
like background melts. 



To investigate superconducting coherence, in Figure ^ we show the inverse dielectric 
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function versus temperature. We see that vanishes at the defect melting transition Tm- 
The diffusing defects above induce a diffusion of vortices, which must move to fill in the 
"hole" left behind by the defect as it moves. The diffusing vortices are then responsible for 
the destruction of superconducting phase coherence. 

To summarize, we have found clear evidence that the introduction of a small concen- 
tration of defects into the fully frustrated system results in a dramatic decrease of the 
superconducting transition temperature from its / = 1/2 value. The fluctuations of the de- 
fect superlattice, on an essentially frozen / = 1/2 like background, result in behavior which 
is in many respects like that of the dilute vortex lattices studied in Section IV. In the present 
case, we find that the defect superlattice melts directly from a pinned sohd into a liquid. 
In the following section we will argue, following the analogy with Section IV, that a more 
dilute defect superlattice would first unpin at a Tc(/) into a floating defect superlattice, 
which would then melt at a higher temperature into a defect liquid. 

The sharp melting transition of the / = 1/2 like background at a temperature T^' 
distinctly higher than the defect melting at is a new phenomenon, with no analogue 
in the dilute small / systems (at small / there is only a smooth crossover remnant of the 
vortex-antivortex unbinding transition of the / = case). From symmetry, one would expect 
that the transition at Tm' is of the Ising type. Undersanding whether the melting transition 
in the pure / = 1/2 case is Ising like or not, has been the subject of much work, with the 
most recent simulations suggesting that it is not;0 if it is not Ising, this is most likely due 
to the long range nature of the vortex interactions. For the / = 5/11 case however, the 
melted defect liquid will serve to screen the interactions of the vortices in the / = 1/2 like 
background, resulting in effectively short ranged interactions. An Ising transition is therefore 
most probable. We are unable to test this prediction, as we are unable to carry out a 
detailed finite size scaling analysis, for the same reasons as discussed in Section IV. However 
the strong screening effect of the defect liquid is evident in the substantial reduction of the 
background melting temperature, Tm'{5/ll) ~ 0.04, as compared to the melting transition 
T^(l/2) ~ 0.13 of the pure / = 1/2 case. 
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B. General case: / = l/2 — l/g 



In this section, we strengthen the analogy between the melting of the defect superlattice 
seen in the previous section, with the melting of the dilute (small /) vortex lattices studied 
in Section IV, in order to discuss the general case of / = 1/2 — Our goal is to establish 
that for a more dilute density of defects, one will have a floating defect solid intermediate 
between a pinned defect solid and a defect liquid. 

As seen for the / = 5/11 case, throughout a temperature range including the defect 
lattice melting transition T^, the / = 1/2 like background vortices remain perfectly ordered 
as at T = 0; domain excitations, which would reduce 5'(q*) from its T = value of unity, 
become important only above Tm- In this case, one can focus solely on excitations which 
are due to the motion of the defects. At T = 0, these defects are seen to sit on the same 
sublattice of the square grid as do the rij = 1 vortices of the / = 1/2 like background 
(see Figure equivalently, one never has two vortices on two nearest neighbor sites. We 
assume that this restriction continues to hold at finite temperatures up to and including Tm, 
i.e. the cost in energy to have two vortices on nearest neighbor sites is so high compared to 
Tm, that such excitations may be ignored. 

With this assumption, we have reduced our problem at / = 1/2 — 1/gto that of a density 
1/g of logarithmically interacting defects, which are restricted to move on only one of the 
two sublattices of the original square grid. This sublattice is itself a square lattice of lattice 
constant \/2aQ. As the sublattice has half the number of sites as in the original grid, our 
problem is thus effectively the same as a dilute density /' = 2/g of vortices on a square 
grid. This mapping would be exact (within our assumptions) except for the fact that the 
interaction potential ^(r) between the defects is still defined with respect to the original 
square grid, and not the sublattice to which the defects are constrained. However, as q gets 
large, and the average spacing between defects becomes much greater than oq, we expect 
that this difference will be a negligible effect. 

The assumption that the defects move only on one sublattice can be checked for the 
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/ = 5/11 case of the previous section. Restricting the defects in real space to a sublattice 
whose unit cell has twice the area of that of the orginial square grid, means that the 1st 
Brillouin Zone of the effective reciprocal lattice is reduced by a factor of one half. Instead 
of the square shaped BZ shown in Figures the effective BZ is now an inscribed diamond 
whose vertices bisect the edges of the squares of Figures pil. The structure function S(k), 



as plotted over the full square shaped BZ of the original square grid, should now just be 
obtained by a periodic repetition of the diamond shaped BZ corresponding to the sublattice. 
Such periodicity is clearly seen in Figures pl| a and |2l|b, for both the pinned defect solid, and 



the melted defect liquid. It is absent in Figure pi|c, where T > T^', and the / = 1/2 like 
background has melted. 

Having checked the validity of our assumption for / = 5/11, we note that that it should 
be even better satisfied for more dilute defect densities / = 1/2 — 1/g, g > 22. As g increases, 
the density of defects decreases, resulting a reduced screening of the interactions between 
the background vortices. The background melting temperature T^' should therefore increase 
and approach its higher / = 1/2 value. At the same time, the defect superlattice unpinning 
temperature Tc should decrease as ~ 1/g, while the defect superlattice melting temperature 
Tm saturates to a lower fixed value. Thus we expect that the window of temperatures in 
which our assumption is valid becomes wider as g increases. 

We can now understand the behavior found in the preceeding section. For / = 5/11 = 
1/2 — 1/22, we have g = 22, and so the defects behave like an effective vortex density of 
/' = 2/g = 1/11. Comparing to our results of Figure |1^ in Section IV, we see that /' is large 
enough that we expect only a direct melting of the pinned defect solid to a defect liquid, 
consistent with our observation in the preceeding section. In order to observe a floating 
defect solid, we will have to consider an /' < 1/30, or an / = 1/2 — 1/g with g > 60. 

To simulate a system with / = 1/2 — 1/60 directly, would require a grid size of at 
least L = 60, with Nc = 1740 vortices. This is beyond our present computational ability 
(/ = 5/11, with L = 22 and Nc = 220, is about the largest system we can manage). However 
our conclusion, that at temperatures low compared to Tm' the defects move in the presence 
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of an effectively frozen / = 1/2 like background, allows us to construct a much more efficient 
algorithm which will be suitable for describing behavior up to and above the defect melting 
transition T^, provided we stay below the background melting transition T^'. We do this 
by fixing the / = 1/2 like background and allowing only the vacancies to move around. This 
significantly reduces the number of degrees of freedom in the simulation, and we shall thus 
be able to treat systems with a much smaller fraction of defects 1/g, than we could by direct 
MC simulation. 

In order to implement the algorithm suggested above, we formally decompose the charge 
at site Ti into two parts 

ni = Si- 6ni, (33) 

where 

= i[l + (-1)-^+^^] = 1(1 + e--'i*) (34) 

is the staggered pattern of the background vortices (q* = (vr/ao, vr/ao)) and 6ni are new 
integer variables representing the defects in the background. Neutrality requires that 
J^i^^i = L'^/q. Substituting Eq.(|55|) for rii in the Hamiltonian (^), we get 

^ = ^ E 5n,V:^Sn, - Y: 5n,Vris, - /) + ^ E(^^ - fW^M ' /)' (35) 

ij ij ij 

where V/j = V'{ri — rj). The first term in the Hamiltonian (|35D gives the interaction between 
defects; the second term represents the interaction of the defects with a one-body potential 

j 

created by the background; the last term is just an additive constant. Substituting Eq.(plD 
for the Si into Eq.(^) above, gives the potential $j in terms of the Fourier components of 
the interaction, Vk, 

<^>. = Iv^^e-^-^* - - /) E ^k, (37) 
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where from Eq.(^ we have V^. = 7r/4. Thus oscillates with the same checkerboard pattern 
of the / = 1/2 like background. Comparing with Eq.(0), we see that the Hamiltonian of 
Eq.(p5D can now be rewritten in the following simple form 

1 TT 

= - ^ 6niV-j6nj - ^Yl ^^i^i + -^o, (38) 

ij i 

where Eq is an additive constant. 

So far, the formulation above is exact. Our approximation that the background is frozen, 
and that defects only move on the sublattice defined by Sj = 1, occurs when we consider 
only the case where Nc sites have the value Srii = 1, and all other sites have Sui = 0. In this 
approximation, the Coulomb gas near full frustration, /=l/2 — is equivalent at low 
temperatures to the dilute Coulomb gas of defects with integer charges 6ni = 0, 1, moving in 
a staggered potential of magnitude 5$ = 7r/4 = 0.7853 .... As this magnitude is about two 
orders of magnitude greater than the relevant excitation energy scale, set by temperature 
Tm, the sites with 6ni = 1 are essentially restricted to the sublattice where Sj = 1; in this 
case they represent the vacancies in the / = 1/2 like background. The case where 6ni = 1 
on the opposite sublattice where Si = 0, represents a (+1, —1) vortex-antivortex excitation, 
which can be ignored on energtic grounds as we had shown in earlier sections. 

To check the consistency of the above procedure, we have redone our simulation of 
/ = 5/11 using the new algorithm based on the Hamiltonian (^). In a fraction of the CPU 
time needed for the original simulation using the full Hamiltonian (^), we have recovered 
our original results for all quantities, at all temperatures up to about T ^ 0.040, where 
fluctuations in the / = 1/2 background become significant. 

Having verified the consistency of the new algorithm in this way, we now proceed to 
simulate systems with more dilute concentrations of defects. In Figure |2^, we display the 
structure function S(k) for the case / = 22/45 = 1/2 — 1/90. As expected from the 
discussion above, we observe a clear signature of the floating solid phase (Figure in 
the temperature range 0.005 < T < 0.008. This range is identical to the range in which we 



found the floating solid phase for /' = 1/45 (see Figure |T^). The low temperature phase is 
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a familiar "pinned defect solid" (Figure [2^a); the high temperature phase is a defect liquid 
with strong four-fold correlations (Figure p4|c). The above scenario is confirmed by a direct 
measurement of e~^(T) and the six-fold orientational correlation ^PqIT) of the defects Srii, 
shown in Figure |2^. Both quantities behave in a way similar to those measured for dilute 
vortex systems, showing a sharp drop in e^^(T) at the depinning transition, and a plateau 
in feiT) in the floating phase. 

To summarize, we conclude that for / = 1/2 — with q > 60, there will be the 
following sequence of transitions. At low temperature there is a pinned superlattice of 
defects of density 1/q, which unpins at Tc(/) ~ 2/g into a floating superlattice of defects. 
This floating lattice melts at — 0.007 into an isotropic defect liquid. Finally, at T^', 
which approaches the value of 0.13 as q increases, the / = 1/2 like background melts via an 
Ising transition, resulting in an isotropic vortex liquid of density /. 

VI. SUMMARY AND CONCLUSIONS 

We have carried out extensive Monte Carlo simulations of the Coulomb gas Hamiltonian 
(^) as a model of a 2D superconducting network in an external transverse magnetic field. One 
of the goals of our work was to systematically study, for the special cases of vortex density 
f = 1/q and / = 1/2 — l/g(g^2)a conjecture put forward by Teitel and Jayaprakash,! that 
for / = p/q the superconducting transition temperature scales approximately as Tc{f) ^ 1/q. 
For the dilute case, f = 1/q, we have found good agreement with this conjecture, provided 
one interprets the superconducting transition temperature to be the vortex lattice unpinning 
temperature Tc{f), where the ground state vortex lattice decouples from the superconducting 
network, and is free to slide transversly to any applied d.c. current, thus producing "flux 
flow" resistance. 

A new result of our work is the realization that above Tc(/), for sufficiently dilute systems, 
a depinned "floating" vortex solid will exist. This floating vortex solid has essentially the 
same properties as a vortex lattice in a uniform superconducting film, and it melts (as 
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q oo) at Tm — 0.007 into an isotropic vortex liquid. While the true onset of finite linear 
d.c. resistivity will be Tc{f), the melting at is presumably accompanied by a sharp rise 
in resistivity. The distinction between Tc and Tm, however, may be difficult to observe 
experimentally, due to the existence of large energy barrierJi^ for the hopping of a vortex 
between neighboring cells of the superconducting network. As discussed in the introduction, 
the discrete nature of the network introduces an effective periodic pinning potential for 
vortices. For a square Josephson array, Lobb et a/.0 have estimated the energy barrier of 
this pinning potential to be Eb — 0.199/27r = 0.0317 (in our energy units0). This is almost 
five times the vortex lattice melting temperature Tm — 0.007! Thus for the square network, 
one is most likely to observe upon cooling only a vortex liquid, in which the vortex mobility 
decreases exponentially as e"^*"/^; the true phase transitions at Tm and will be masked by 
the extremely slow relaxation over the energy barriers Ef, at these low temperatures. Such 
behavior has in fact been reporte in experimental studies of square Josephson arrays, 
where for small / near / = 0, only exponentially decreasing resistive tails are observed 
at low temperature; no evidence for the melting or depinning transitions at Tm and T^ 
has been found. For the triangular Josephson array however (a case we have not explicitly 
studied here), the energy barrier is estimatedB to be Ej, ~ 0.0427/27r = 0.0068 (in our energy 
units). This is comparable to Tm, and so there might be some slight chance of experimentally 
observing the melting transition. Recent experimental studies of this systemil at small / 
have found surprising dynamical behavior, indicating anomalously slow diffusion of vortices. 
However the temperature of these experiments, T ~ 0.5Tc(/ = 0), appears to be too high for 
these results to be explained by any of the melting or depinning effects we have found here. 
For a honeycomb Josephson array (vortices on a triangular grid) we estimate the highest 
barrier,^ E^ ~ 0.751/27r = 0.119. 

For dense systems close to full frustration, / = 1/2 — 1/g, we have argued that, at 
temperatures low compared to the melting temperature of the pure / = 1/2 system, the 
physics is dominated by defects moving on top of a quenched / = 1/2 like vortex background. 
We have shown how this system of defects can be mapped onto the dilute Coulomb gas of 
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vortices with density /' = 2/q. The resulting behavior is then obtained from knowledge of 
this dilute limit. The TJ conjecture again holds, with Tc(/ = (g — 2)/2g) ~ 2/g marking the 
transition from a pinned defect superlattice to a floating defect superlattice, in which true 
d.c. superconductivity is lost. At a higher — 0.007, the floating defect solid melts into an 
isotropic defect liquid. At yet a higher T^', the / = 1/2 like background disorders. As in the 
/ = case, the transitions at and Tm may be difficult to observe experimentally, due to the 
energy barrier for a defect to hop between nearest neighbor sites of the relevant sublattice. 
For a square lattice. Dang and Gyorffyii have estimated this barrier to be Ei, ~ 0.368/27r 
(in our energy units), even larger than that found for / = 0. Experimental studieJllii of 
square Josephson arrays for / near / = 1/2 have again found only exponential resistive tails 
as the temperature decreases. 

In contrast to the transitions at and T^, we expect that the sharp disordering transi- 
tion of the / = 1/2 like background at T^' should be experimentally observable. This follows 
since for / = 1/2 — 1/g, we expect that as g oo, T^' Tm(l/2) ~ 0.13, well above the 
energy scale of the barriers. This transition would presumably manifest itself as a singular 
increase in the linear resistivity at T^/ (from an already finite value). The phase boundary 
Tm'if) near / = 1/2 is presumably a smooth funtion of the defect density, 1/2 — f, however 
we are unable to estimate it due to our inability to simulate sufficiently large systems. 

Our mapping between a dilute density of defects near / = 1/2, and a dilute vortex 
density near / = 0, may be extended to the more general case. Using the same arguements 
as in Section V, we would expect that the system with / = 1/2 —p/q, with p/q sufficiently 
small, should have the same low temperature behavior as the density /' = 2p/q. For p/q 
sufficiently small, we would expect that the dilute vortex lattice of density /' = 2p/q behaves 
qualitatively like those of density 1/q studied here, i.e. there is first a depinning transition 
at Tc(/) which decreases asp/g decreases (whether it vanishes asp/g or as 1/g remains to be 
investigated) followed by a melting transition at Tm ~ 0.007. Thus for any rational fraction 
sufficiently close to either / = or / = 1/2, we would expect behavior similar to the cases 
explicitly studied here. 
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We thus see that the TJ conjecture appears to hold, according to the following scenario. 
Consider a vortex density / close to some simple fraction /q = Po/qo, f = /q — with go 
q. The ground state is one which is almost everywhere like that of /o, except for a pinned 
periodic superlattice of defects of density 1/q. If g is sufficiently large, this superlattice will 
unpin into a floating defect solid at Tc{f) ^ 1/q. Defects which are free to move lead to flux 
flow resistance, and destroy the superconducting phase coherence of the system. Thus an 
arbitrarily small concentration of defects added on top of the /g like background (i.e. for / 
arbitrarily close to /o) dramatically decreases the superconducting transition temperature, 
when compared to pure /o system. We have explicitly tested this scenario for the cases /o = 
and /o = 1/2. We speculate that this behavior will be characteristic of systems near any 
simple fraction /q = Po/qo- We further speculate that his behavior may be characteristic 
for any sufficiently small rational fraction of defects away from a simple fraction /q, i.e. 
/ = Po/qo -p/q with go < q- 

A second goal of our work was to study in detail the melting transition of the 2D vortex 
lattice. This problem has been addressed previously only in the context of uniform super- 
conducting films. Here we have addressed this issue in the context of a superconducting 
network. We believe, however, that our results for the dilute case we have studied in Section 
III are representative of the continuum limit, as treated within the London approximation. 
Melting within this London approximation has been treated theoretically by Huberman and 
Doniach ,§ and FisherS, who applied the KTNHY theory of defect mediated melting in 2D. 
This theory predicts a second order melting transition at a well below the Ginzburg- 
Landau mean field transition temperature Tmf, as well as an intermediate hexatic liquid 
phase. We have carried out the first detailed finite size scaling analysis to check this KTNHY 
theory as applied to vortex lattice melting. We find a value — 0.007 ± 0.0005 in good 
agreement with the value estimated by Fisher. We also find that the vortex lattice shear 
modulus jumps discontinuously to zero at T^, with a value close to the KTNHY prediction. 
However we find that the melting transition is weakly first order, and we find no evidence 
for a hexatic phase. Our value for T^, and our conclusion concerning the order of the melt- 
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ing transition, are in agreement with earlier simulationsEl of the continuum one component 
plasma model of Eq.(|l^). 

This problem of 2D vortex lattice melting has been the focus of much renewed work 
recently, due to its potential connection with behavior in anisotropic high temperature su- 
perconductors. The very existence of a vortex lattice at any finite temperature has been 
challenged by Moore,i who argued that fluctuations in the phase of the order parameter 
ip{r), due to shear excitations of the vortex lattice, will cause the order parameter corre- 
lation function {ip*{r)ip{0)) to decay exponentially at any finite temperature. From such 
decay, Moore argued first for the absence of a superconducting state, and then concluded 
as a result of this absence of superconductivity, that the vortex lattice should not exist. 
Support for this scenario is suggested by high temperature perturbative expansions, which 
also find no evidence for freezing into a vortex lattice, even when evaluated to high order.@ 

Recently, simulations have been carried out to address this question. In contrast to our 
work in the London approximation, these works have been carried out in the so called lowest 
Landau level (LLL) approximation, which is based on the Ginzburg-Landau (GL) free energy 
TicL of Eq.(|TT]). In this approximation, the complex order parameter ip{r) is expanded in 
terms of the lowest degenerate eigenstates of the Gaussian part of the Ginsburg-Landau 
free energy, and the coefficients of this expansion (or alternatively the complex positions 
of the vortices) are used as fluctuating variables in a Monte Carlo simulation with TicL as 
the Hamiltonian. Using such simulations, and modeling a 2D system by the surface of a 
sphere, O'Neill and Moorell^ failed to find evidence for a vortex lattice. Other simulations in 
a 2D plane however reported clear evidence for the melting of a vortex lattice at a finite 
temperature. Hu and MacDonaldlil and Kato and NagaosaH find that this melting transition 
is weakly first order, in agreement with our London result. Sasik and Stroud^ similarly find 
a first order transition; they further compute the vortex lattice shear modulus fi and find 
behavior at Tm in agreement with our result. Most recently, Herbut and Tesanovic^ have 
developed a density functional theory of the vortex lattice melting transition, based on the 
LLL formalism. They again find results consistent with the above, for the order of the 
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transition, and the shear modulus /i. 



Thus, with the exception of Ref. ^ results from the London and LLL approximations 
seem to be in agreement. This is as one might expect from the principle of universality in 
phase transitions. Although the London approximation at the "microscopic" level ignores 
fluctuations in the amplitude of ip{r) (such as are included in the LLL formalism), upon 
coarse graining the London model, phase fluctuations at the microscopic length scale will 
generate amplitude fluctuations on the coarse grained length scale. On this coarse grained 
scale, the system will be described by some effective Ginzburg-Landau free energy, complete 
with amplitude fluctuations, although higher order terms in ip beyond those given in Eq.(|TTD 
may be present. In contrast to the London approximation, the Ginzburg-Landau free energy 



of Eq.(pJ]) includes amplitude fluctuatons at the "microscopic" scale, and it is thus often 
viewed as a more fundamental description. However, the GL form of Eq. (0) represents only 
the lowest terms in an expansion of the free energy in powers of and hence is valid only 
near the mean field transition Tmf where is small. Since vortex lattice melting occurs at 
<< Tmf, higher order terms in ip may well be important for a quantitative description. 
These higher order terms, however, are presumably irrelevent in determining the critical 
behavior, thus leading to agreement between the London and LLL simulations. 
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TABLES 



T 








0.00475 


0.188 ±0.008 


0.704 ±0.055 


0.571 ±0.007 


0.00500 


0.207 ±0.007 


0.806 ±0.032 


0.529 ±0.005 


0.00525 


0.211 ±0.007 


0.852 ±0.028 


0.504 ±0.004 


0.00550 


0.248 ±0.005 


0.998 ±0.019 


0.476 ±0.003 


0.00575 


0.255 ±0.008 


0.999 ±0.029 


0.458 ±0.003 


0.00600 


0.296 ±0.006 


1.065 ±0.028 


0.426 ±0.007 


0.00625 


0.319 ±0.010 


1.191 ±0.016 


0.403 ±0.004 


0.00650 


0.4 ±0.16 


1.4 ±0.22 


0.33 ±0.030 


0.00675 


1.4 ±0.31 


2.0 ±0.31 


0.20 ±0.041 


0.00750 


3.4 ±0.37 


3.4 ±0.44 


0.03 ±0.046 


0.01100 


2.8 ±0.23 


2.9 ±0.30 


-0.01 ±0.032 


0.01500 


2.2 ±0.12 


2.1 ±0.22 


0.00 ±0.020 



TABLE I. Temperature dependence of the exponents r/Oi {T) and 77G2 {T) for / = 1/49 on the 
triangular grid, as obtained from finite size scaling. Also displayed are the limiting values yjg" of 
the orientational correlation ^q{T) for L — > 00. 
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T 




L = 63 


L = 77 


L = 91 


FSS 


0.00475 


0.164±0.026 


0.161±0.025 


0.195±0.018 


0.188 ±0.008 


0.00500 


0.216±0.012 


0.221±0.009 


0.219±0.007 


0.207 ±0.007 


0.00525 


0.249±0.009 


0.235±0.009 


0.247±0.006 


0.211 ±0.007 


0.00550 


0.282±0.007 


0.272±0.006 


0.275±0.007 


0.248 ±0.005 


0.00575 


0.298±0.008 


0.292±0.009 


0.290±0.006 


0.255 ±0.008 


0.00600 


0.326±0.008 


0.326±0.007 


0.329±0.007 


0.296 ±0.006 


0.00625 


0.351±0.014 


0.350±0.017 


0.352±0.016 


0.319 ±0.010 


0.00650 


0.677±0.342 


0.473±0.223 


0.568±0.131 


0.4 ±0.16 


0.00675 


1.755±0.789 


1.678±0.453 


2.458±0.911 


1.4 ±0.31 



TABLE II. Comparison of the exponents r/d (T) obtained using two different methods. 
Coulumns 2-4 show results from fitting of S{G) to Eg. (p4[) , for system sizes L = 63, 77, 91. Coul- 
umn 5, labelled FSS, restates the results from the finite size scaling analysis. All exponents are for 
density / = 1/49. 
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T 




0.0040 


0.0024±0.001 


0.0045 


0.111 ±0.016 


0.0050 


0.198 ±0.012 


0.0055 


0.224 ±0.009 


0.0060 


0.270 ±0.011 


0.0065 


0.33 ±0.04 


0.0070 


0.49 ±0.10 



TABLE III. Temperature dependence of the exponents r/Ci (T) of the floating vortex lattice 
on the square grid as obtained by fitting the height of peaks in the structure function, for / = 1/60 
and L = 60. 
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FIGURES 

FIG. 1. Phase diagram of the sufficiently dilute system, as found by our Monte Carlo calcu- 
lation. 

FIG. 2. Structure function S'(k) in the first Brillouin zone, (upper portion of the figure), and 
the profile of the peak heights along the vertical symmetry axis ky (lower portion); for / = 1/49 
and Nc = 63, and three different temperatures T: (a) T = 0.003, just below Tc, in the "pinned 
solid" phase, (b) T = 0.0065, just below Tm, in the "floating solid" phase, (c) T = 0.0075, just 
above Tm, in the liquid. Intensities in the density plots are plotted nonlinearly to enhance features. 

FIG. 3. Inverse dielectric function e^^(T) and orientational order correlation ^&{T) versus T 
for / = 1/49 and = 169. Solid and dashed lines are guides to the eye only. 

FIG. 4. The dependence of the depinning and melting temperatures, Tc and T^, on vortex 
density /. Errors are estimated from the width in the apparent drop in e^^(r) and ipQ{T). Solid 
and dashed lines are guides to the eye only. 

FIG. 5. Finite size scaling of S'(Gi)/L^ (note the log-log scale) for / = 1/49. Solid and dashed 
lines are fits to Eq. (p2|b). 

FIG. 6. Heights of the peaks ^(G) versus |G| for / = 1/49 and L = 63. Dashed lines represent 
the best fit to Eq. (p4|) , and are used to extract the exponent 7/Gi(T). 

FIG. 7. Finite size scaling of (feiT) for / = 1/49. Solid and dashed lines are fits to Eq. (p7|a). 

FIG. 8. T/r}G^{T) (proportional to the shear modulus /i) and orientational correlation 93g° 
versus T, as extracted from finite size scaling. The intersection of T/rja^ (T) with the dashed line 
3r determines the KTNHY upper bound on the melting transition T/r]G-^(T) > 3r. 
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FIG. 9. Comparison of Vcli"^) ^nd -l/K\n[(p^{T)] as a test of Eq.(^). Below the two 
quantities coincide, as expected for a 2D harmonic lattice. Dashed line determines the KTNHY 
upper bound t^q^ (T) > 3. 

FIG. 10. Translational and orientational correlation lengths $.+(T) and (,q{T) versus T for 
/ = 1/49, as extracted from finite size scaling. Both correlation lengths sharply increase at the 
melting temperature Tm ~ 0.0070. 

FIG. 11. Free energy distribution F[E) versus E, at melting T^, for / = 1/49 and several 
system sizes L. The growth in the energy barrier AF with increasing L (see inset) indicates a first 
order transition. Curves for different L are offset from each other by a constant, for the sake of 
clarity. 

FIG. 12. Free energy distribution F{E) versus E, at the depinning transition Tc, for / = 1/49 
and several system sizes L. The growth in energy barrier AF with increasing L (see inset) indicates 
a first order transition. Curves for different L are offset from each other by a constant, for the sake 
of clarity. The abrupt ending of the distributions at the low energy side of the graph is because 
the lower minimum represents the ground state energy. 

FIG. 13. Dependence of Tc and Tm on vortex density / for the dilute system on a square grid. 
Dashed and solid lines are guides to the eye only. 

FIG. 14. Two types of ground state configurations for a dilute system on the square grid: (a) 
nearly triangular vortex lattice / = 1/60; (b) nearly square vortex lattice / = 1/51. Solid squares 
denote positions of vortices. 

FIG. 15. (a) Inverse dielectric constant e~^(T) and (b) orientational correlations ^peiT) and 
(pi{T) versus T, for the system with nearly triangular ground state with / = 1/60 and L = 60. 
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FIG. 16. Melting of a nearly triangular vortex lattice on the square grid. Intensity plots of S'(k) 
for / = 1/60, L = 60 and several temperatures: (a) T = 0.003 in the pinned solid; (b) T = 0.006 
floating solid; (c) T = 0.009 in the liquid. 

FIG. 17. Heights of the peaks >S'(G) versus |G| for / = 1/60 and L = 60. Dashed lines 
represent the best fit to Eq. (p4[) , and are used to extract the exponent riQ,^{T). 

FIG. 18. (a) Inverse dielectric constant e~^(T) and (b) orientational correlations ^eiT) and 
9?4(r) versus T, for the system with nearly square ground state with / = 1/51 and L = 51. 

FIG. 19. Melting of a nearly square vortex lattice on the square grid. Intensity plots of ^(k) 
for / = 1/51, L = 51 and several temperatures: (a) T = 0.003 in the pinned solid; (b) T = 0.0045 
in the floating solid; (c) T = 0.006 in the liquid. 

FIG. 20. Ground state for / = 5/11 on a 22 x 22 unit cell. Solid squares represent vortex 
positions. Crosses (+) indicate vacancies (defects) in the otherwise perfect checkerboard pattern 
of / = 1/2. 

FIG. 21. Melting of / = 5/11 for L = 22. S{k) is shown for: (a) T = 0.010 in the pinned 
defect solid; (b) T = 0.018 in the defect liquid; (c) T = 0.055 in the completely disordered high 
temperature phase. 

FIG. 22. Peak heights S{q)* with q* = (vr/oo, vr/ao), and S{Gi) where Gi = (27r/L)(l,5) is 
the shortest reciprocal lattice vector of the defect superlattice, plotted versus T for / = 5/11 and 
L = 22. 

FIG. 23. Inverse dielectric constant e"^(T) versus T for / = 5/11 and L = 22. 
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FIG. 24. Melting of the defect superlattice for / = 22/45 and L = 90. Intensity plot of ^(k) 
for: (a) T = 0.0040 in the pinned defect solid; (b) T = 0.0065 in the defect floating solid; (c) 
T = 0.0085 in the defect liquid. Intensities at (±7r/ao, ±7r/ao) and (0,0), that arise from the fixed 
staggered background, have been substracted for the sake of clarity. 

FIG. 25. Inverse dielectric function e~^{T) and orientational correlation ipe{T) versus T, for 
/ = 22/45 and L = 90. 
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